################################################################################################%
# This script takes the bootstrapped SV data under the rules specified in the SV Welfare  
# Simulations Script and then produces graphs for some welfare measures as in 
# Figures 4 and 5 and Table 2 of the Columbia Elections Paper.
# 
# All graphs are exported to the ./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/ .
# 
# Author: Luis S.
# Last Modified: 8.25.17
################################################################################################%

library(dplyr)
library(ineq)
library(ggplot2)
library(stargazer)

# setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results/")

SV_EducBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_EducBSPref.csv", header = T, stringsAsFactors = F)
SV_ImmBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_ImmBSPref.csv", header = T, stringsAsFactors = F)
SV_BondBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_BondBSPref.csv", header = T, stringsAsFactors = F)
SV_TeachBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_TeachBSPref.csv", header = T, stringsAsFactors = F)

SV_PrefSummary_raw <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_PrefPointSummary.csv", header = T, stringsAsFactors = F)

SV_RuleAResults <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_RuleAResults.csv", header = T, stringsAsFactors = F)
SV_RuleBResults <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_RuleBResults.csv", header = T, stringsAsFactors = F)
SV_RuleCResults <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_RuleCResults.csv", header = T, stringsAsFactors = F)
SV_RuleDResults <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_RuleDResults.csv", header = T, stringsAsFactors = F)
SV_RuleEResults <- read.csv("./SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_RuleEResults.csv", header = T, stringsAsFactors = F)


SV_RuleAResults <- mutate(SV_RuleAResults, Rule = "A")
SV_RuleBResults <- mutate(SV_RuleBResults, Rule = "B")
SV_RuleCResults <- mutate(SV_RuleCResults, Rule = "C")
SV_RuleDResults <- mutate(SV_RuleDResults, Rule = "D")
SV_RuleEResults <- mutate(SV_RuleEResults, Rule = "E")

SV_AllRules <- bind_rows(SV_RuleAResults, SV_RuleBResults, SV_RuleCResults, SV_RuleDResults, SV_RuleEResults)

####################################
###### Determine Minority Win ######
###################################%


SV_AllRules <- SV_AllRules %>% 
  mutate(SV1_MinorityWon = ifelse(SV1_Winner_BilingualEduc == "Minority" | SV1_Winner_Immigration == "Minority" | 
                                    SV1_Winner_PublicBonds == "Minority" | SV1_Winner_TeacherTenure == "Minority",T, F),
         SV2_MinorityWon = ifelse(SV2_Winner_BilingualEduc == "Minority" | SV2_Winner_Immigration == "Minority" | 
                                    SV2_Winner_PublicBonds == "Minority" | SV2_Winner_TeacherTenure == "Minority",T, F),
         MinorityWinIsEff = ifelse(NSV_Result_BilingualEduc != Eff_Result_BilingualEduc | NSV_Result_Immigration != Eff_Result_Immigration |
                                     NSV_Result_PublicBonds != Eff_Result_PublicBonds | NSV_Result_TeacherTenure != Eff_Result_TeacherTenure, T, F)) %>%
  mutate(SV1IsFullyEff = ifelse(SV1_Result_BilingualEduc != Eff_Result_BilingualEduc | SV1_Result_Immigration != Eff_Result_Immigration |
                                  SV1_Result_PublicBonds != Eff_Result_PublicBonds | SV1_Result_TeacherTenure != Eff_Result_TeacherTenure, F, T),
         SV2IsFullyEff = ifelse(SV2_Result_BilingualEduc != Eff_Result_BilingualEduc | SV2_Result_Immigration != Eff_Result_Immigration |
                                  SV2_Result_PublicBonds != Eff_Result_PublicBonds | SV2_Result_TeacherTenure != Eff_Result_TeacherTenure, F, T),
         SV1IsEff_BE = ifelse(SV1_Result_BilingualEduc == Eff_Result_BilingualEduc,T,F),
         SV1IsEff_IM = ifelse(SV1_Result_Immigration == Eff_Result_Immigration,T,F),
         SV1IsEff_PB = ifelse(SV1_Result_PublicBonds == Eff_Result_PublicBonds,T,F),
         SV1IsEff_TT = ifelse(SV1_Result_TeacherTenure == Eff_Result_TeacherTenure,T,F),
         SV2IsEff_BE = ifelse(SV2_Result_BilingualEduc == Eff_Result_BilingualEduc,T,F),
         SV2IsEff_IM = ifelse(SV2_Result_Immigration == Eff_Result_Immigration,T,F),
         SV2IsEff_PB = ifelse(SV2_Result_PublicBonds == Eff_Result_PublicBonds,T,F),
         SV2IsEff_TT = ifelse(SV2_Result_TeacherTenure == Eff_Result_TeacherTenure,T,F)) %>%
  mutate(SV1_MinorityWon_BE = ifelse(SV1_Winner_BilingualEduc == "Minority",T,F), 
         SV1_MinorityWon_IM = ifelse(SV1_Winner_Immigration == "Minority",T,F),
         SV1_MinorityWon_PB = ifelse(SV1_Winner_PublicBonds == "Minority",T,F),
         SV1_MinorityWon_TT = ifelse(SV1_Winner_TeacherTenure == "Minority",T, F),
         SV2_MinorityWon_BE = ifelse(SV2_Winner_BilingualEduc == "Minority",T,F), 
         SV2_MinorityWon_IM = ifelse(SV2_Winner_Immigration == "Minority",T,F),
         SV2_MinorityWon_PB = ifelse(SV2_Winner_PublicBonds == "Minority",T,F),
         SV2_MinorityWon_TT = ifelse(SV2_Winner_TeacherTenure == "Minority",T, F),
         MinorityWinIsEff_BE = ifelse(NSV_Result_BilingualEduc != Eff_Result_BilingualEduc,T,F), 
         MinorityWinIsEff_IM = ifelse(NSV_Result_Immigration != Eff_Result_Immigration,T,F), 
         MinorityWinIsEff_PB = ifelse(NSV_Result_PublicBonds != Eff_Result_PublicBonds,T,F), 
         MinorityWinIsEff_TT = ifelse(NSV_Result_TeacherTenure != Eff_Result_TeacherTenure, T, F)) %>%
  mutate(EffMinorityWinWithSV1_BE = ifelse(SV1_MinorityWon_BE & MinorityWinIsEff_BE, T, F),
         EffMinorityWinWithSV1_IM = ifelse(SV1_MinorityWon_IM & MinorityWinIsEff_IM, T, F),
         EffMinorityWinWithSV1_PB = ifelse(SV1_MinorityWon_PB & MinorityWinIsEff_PB, T, F),
         EffMinorityWinWithSV1_TT = ifelse(SV1_MinorityWon_TT & MinorityWinIsEff_TT, T, F),
         EffMinorityWinWithSV2_BE = ifelse(SV2_MinorityWon_BE & MinorityWinIsEff_BE, T, F),
         EffMinorityWinWithSV2_IM = ifelse(SV2_MinorityWon_IM & MinorityWinIsEff_IM, T, F),
         EffMinorityWinWithSV2_PB = ifelse(SV2_MinorityWon_PB & MinorityWinIsEff_PB, T, F),
         EffMinorityWinWithSV2_TT = ifelse(SV2_MinorityWon_TT & MinorityWinIsEff_TT, T, F)) %>%
  mutate(EffMinorityWinWithSV1 = ifelse(SV1_MinorityWon & MinorityWinIsEff, T, F), EffMinorityWinWithSV2 = ifelse(SV2_MinorityWon & MinorityWinIsEff, T, F))


#---- realized share of surplus (RSS)

# RSS_SV = (W_SV - R)/(W_Eff - R) where R = sum over all i,k preference intensities divided by 2 (to serve as a welfare floor)

welfFloor <- rowSums(SV_PrefSummary_raw %>% select(-matches("subjects")) %>% select(-matches("Round")), na.rm = TRUE)/2
welfFloor <- data.frame(Round = 1:10000, WelfFloor = welfFloor)

SV_AllRules <- full_join(SV_AllRules, welfFloor, by = "Round") 
SV_AllRules <- SV_AllRules %>% mutate(RSS_Maj = (AggregateWelfare_noBV - WelfFloor)/(AggregateWelfare_Eff - WelfFloor),
                                      RSS_SV1 = (AggregateWelfare_SV1 - WelfFloor)/(AggregateWelfare_Eff - WelfFloor),
                                      RSS_SV2 = (AggregateWelfare_SV2 - WelfFloor)/(AggregateWelfare_Eff - WelfFloor))

#---

SV_RuleAResults <- filter(SV_AllRules, Rule == "A")
SV_RuleBResults <- filter(SV_AllRules, Rule == "B")
SV_RuleCResults <- filter(SV_AllRules, Rule == "C")
SV_RuleDResults <- filter(SV_AllRules, Rule == "D")
SV_RuleEResults <- filter(SV_AllRules, Rule == "E")

bootstrapIterations <- nrow(SV_RuleAResults)


# Print Average Votes and Points InFavor and Against across simulations for each issue
stargazer(as.data.frame(as.matrix(colMeans(SV_PrefSummary_raw[,-1]))) %>% tibble::rownames_to_column("Category") %>% rename(Mean = V1), 
          title = "Average of Votes and Points Per Issue (By Side)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Average of Votes and Points Per Issue (By Side).txt")



######################################
##### Realized Share of Surplus ######
#####################################%



getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}


RSS_summary <- SV_AllRules %>% group_by(Rule) %>% 
  summarise(AvgRSS_Maj = mean(RSS_Maj), Maj_lower95CI = getNthEmpPercentile(RSS_Maj,2.5), Maj_upper95CI = getNthEmpPercentile(RSS_Maj,97.5),
            AvgRSS_SV1 = mean(RSS_SV1), SV1_lower95CI = getNthEmpPercentile(RSS_SV1,2.5), SV1_upper95CI = getNthEmpPercentile(RSS_SV1,97.5), 
            AvgRSS_SV2 = mean(RSS_SV2), SV2_lower95CI = getNthEmpPercentile(RSS_SV2,2.5), SV2_upper95CI = getNthEmpPercentile(RSS_SV2,97.5))


RSS_Test <- SV_AllRules %>% 
  mutate(Diff_SV1minusMaj = RSS_SV1 - RSS_Maj,
         Diff_SV2minusMaj = RSS_SV2 - RSS_Maj) %>%
  group_by(Rule) %>% 
  summarise(AvgDiff1_SV1minusMaj = mean(Diff_SV1minusMaj), Diff1_lower95CI = getNthEmpPercentile(Diff_SV1minusMaj,2.5), Diff1_upper95CI = getNthEmpPercentile(Diff_SV1minusMaj,97.5), 
            AvgDiff2_SV2minusMaj = mean(Diff_SV2minusMaj), Diff2_lower95CI = getNthEmpPercentile(Diff_SV2minusMaj,2.5), Diff2_upper95CI = getNthEmpPercentile(Diff_SV2minusMaj,97.5))




sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Realized Share of Surplus (By Rule).txt")
stargazer(as.data.frame(RSS_summary), title = "Realized Share of Surplus (RSS)", type = "text", digits = 4, summary = FALSE, rownames = FALSE)
stargazer(as.data.frame(RSS_Test), title = "Tests of Differences in RSS", type = "text", digits = 4, summary = FALSE, rownames = FALSE)
sink()




##################################################
####### Debug Table (per issue & per rule) #######
#################################################%

# 1. Number of simulations won by the InFavor side under majority
# 2. Number of sims won by the InFavor side under efficiency
# 3. Number of sims in which majority is efficient.
# 4. Number of sims won by the InFavor side under SV
# 5. Number of sims in which SV is efficient.
# 6. Number of sims in which SV is equivalent to majority.
# 7. Number of sims in which SV differs from maj and the InFavor side wins under SV.
# 8. Number of efficient minority victories under SV. 

ruleList <- c("A", "B", "C", "D", "E")
propList <- c("Immigration", "BilingualEduc", "PublicBonds", "TeacherTenure")
abbrevList <- c("IM", "BE", "PB", "TT")


finalResults <- data.frame()
# currentRule <- "C"


for(currentRule in ruleList)
{
  for(index in 1:4)
  {
    # index <- 3
    currentProp <- propList[index] 
    currentAbbrev <- abbrevList[index] 
    
    tempResults <- data.frame(Rule = currentRule, Proposal = currentProp, stringsAsFactors = F)
    
    # Number of simulations won by the InFavor side under majority
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("NSV_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("NSV_Result_",currentProp," == 'InFavor'", sep = "")) 
    tempResults$nSims_InF_Maj <- temp$Count
    
    # Number of sims won by the InFavor side under efficiency
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("Eff_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("Eff_Result_",currentProp," == 'InFavor'", sep = "")) 
    tempResults$nSims_InF_Eff <- temp$Count
    
    # Number of sims won by the InFavor side under, say, SV, rule A.
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("SV1_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("SV1_Result_",currentProp," == 'InFavor'", sep = "")) 
    tempResults$nSims_InF_SV1 <- temp$Count
    
    # Number of sims in which majority is efficient.
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("MinorityWinIsEff_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MinorityWinIsEff_",currentAbbrev," == 'FALSE'", sep = "")) 
    tempResults$nSims_MajIsEff <- temp$Count
    
    # Number of sims in which SV is efficient.
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("SV1IsEff_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("SV1IsEff_",currentAbbrev," == 'TRUE'", sep = "")) 
    tempResults$nSims_SV1IsEff <- temp$Count
    
    # Number of sims in which SV is equivalent to majority.
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% 
      mutate(MajEqualsSV_BE = ifelse(NSV_Result_BilingualEduc == SV1_Result_BilingualEduc,T,F),
             MajEqualsSV_IM = ifelse(NSV_Result_Immigration == SV1_Result_Immigration,T,F),
             MajEqualsSV_PB = ifelse(NSV_Result_PublicBonds == SV1_Result_PublicBonds,T,F),
             MajEqualsSV_TT = ifelse(NSV_Result_TeacherTenure == SV1_Result_TeacherTenure,T,F)) %>%
      group_by_(paste("MajEqualsSV_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MajEqualsSV_",currentAbbrev," == 'TRUE'", sep = "")) 
    tempResults$nSims_MajEqualsSV1 <- temp$Count
    
    # Number of sims in which SV differs from maj and the InFavor side wins under SV rule A.
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% 
      mutate(MajEqualsSV_BE = ifelse(NSV_Result_BilingualEduc == SV1_Result_BilingualEduc,T,F),
             MajEqualsSV_IM = ifelse(NSV_Result_Immigration == SV1_Result_Immigration,T,F),
             MajEqualsSV_PB = ifelse(NSV_Result_PublicBonds == SV1_Result_PublicBonds,T,F),
             MajEqualsSV_TT = ifelse(NSV_Result_TeacherTenure == SV1_Result_TeacherTenure,T,F)) %>%
      group_by_(paste("MajEqualsSV_",currentAbbrev, sep = ""), paste("SV1_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MajEqualsSV_",currentAbbrev," == 'FALSE'", sep = ""), 
                                         paste("SV1_Result_",currentProp," == 'InFavor'", sep = "")) 
    tempResults$nSims_MajDiffSV1_SVInF <- temp$Count
    
    
    # Number of efficient minority victories under SV rule A. 
    temp <- SV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("EffMinorityWinWithSV1_",currentAbbrev, sep = "")) %>%
      summarise(Count = n())%>% filter_(paste("EffMinorityWinWithSV1_",currentAbbrev," == 'TRUE'", sep = "")) 
    tempResults$nSims_EffMinWinSV1 <- temp$Count
    
    
    # Bind results
    finalResults <- bind_rows(finalResults, tempResults)
  }
}

write.csv(finalResults, "SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Summary of Sims (per issue & per rule).csv", row.names = F)



##############################################################
####### Prop of Victories In Favor of Issue (per rule) #######
#############################################################%

# NSV

temp1 <- SV_AllRules %>% group_by(Rule, NSV_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NSV_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- SV_AllRules %>% group_by(Rule, NSV_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NSV_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- SV_AllRules %>% group_by(Rule, NSV_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NSV_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- SV_AllRules %>% group_by(Rule, NSV_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NSV_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

SimpleMajSummary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
SimpleMajSummary <- SimpleMajSummary %>% mutate(Vote = "Simple Majority")


# SV1

temp1 <- SV_AllRules %>% group_by(Rule, SV1_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV1_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- SV_AllRules %>% group_by(Rule, SV1_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV1_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- SV_AllRules %>% group_by(Rule, SV1_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV1_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- SV_AllRules %>% group_by(Rule, SV1_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV1_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

Vote1Summary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
Vote1Summary <- Vote1Summary %>% mutate(Vote = "SV1")



# SV2

temp1 <- SV_AllRules %>% group_by(Rule, SV2_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV2_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- SV_AllRules %>% group_by(Rule, SV2_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV2_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- SV_AllRules %>% group_by(Rule, SV2_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV2_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- SV_AllRules %>% group_by(Rule, SV2_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(SV2_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

Vote2Summary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
Vote2Summary <- Vote2Summary %>% mutate(Vote = "SV2")


# Combine and Print

sink("SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Prop of Victories In Favor of Each Proposal (per Rule & Vote).txt")
stargazer(as.data.frame(SimpleMajSummary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(Vote1Summary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(Vote2Summary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()



##############################################
####### Freq Histogram of Minority Win #######
#############################################%

#### Frequency of elections in which at least one minority victory occurs per rule (Fig 4 in Elections paper). ####%

# SV1
SV1_AllRules <- group_by(SV_AllRules, Rule, SV1_MinorityWon) %>%
  summarize(Freq = round(n()/bootstrapIterations, digits = 3)*100) %>%
  filter(SV1_MinorityWon == TRUE)

welfarePlot <- ggplot(data = SV1_AllRules, aes(x=Rule, y = Freq)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Frequency of Minority Victories") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV1 (No Recal) - Bootstrapped Freq. of Minority Victories")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV1 (No Recal) - Freq. Minority Victories.png", height = 5.5, width = 5.5)


# SV2
SV2_AllRules <- group_by(SV_AllRules, Rule, SV2_MinorityWon) %>%
  summarize(Freq = round(n()/bootstrapIterations, digits = 3)*100) %>%
  filter(SV2_MinorityWon == TRUE)

welfarePlot <- ggplot(data = SV2_AllRules, aes(x=Rule, y = Freq)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Frequency of Minority Victories") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV2 (No Recal) - Bootstrapped Freq. of Minority Victories")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV2 (No Recal) - Freq. Minority Victories.png", height = 5.5, width = 5.5)


###################################################
####### Freq Table of Minority Win (byRule) #######
##################################################%

#### Frequency Table of elections in which at least one minority victory occurs per rule (Table 2 in Elections paper). ####%

# SV1
Summary_Immigration <- group_by(SV_AllRules, Rule, SV1_Winner_Immigration) %>%
  summarize(SV1_Winner_Immigration_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV1_Winner_Immigration == "Majority") %>% mutate(SV1_Winner_Immigration_Freq = 100 - SV1_Winner_Immigration_Freq)

Summary_PublicBonds <- group_by(SV_AllRules, Rule, SV1_Winner_PublicBonds) %>%
  summarize(SV1_Winner_PublicBonds_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV1_Winner_PublicBonds == "Majority") %>% mutate(SV1_Winner_PublicBonds_Freq = 100 - SV1_Winner_PublicBonds_Freq)

Summary_TeacherTenure <- group_by(SV_AllRules, Rule, SV1_Winner_TeacherTenure) %>%
  summarize(SV1_Winner_TeacherTenure_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV1_Winner_TeacherTenure == "Majority") %>% mutate(SV1_Winner_TeacherTenure_Freq = 100 - SV1_Winner_TeacherTenure_Freq)

Summary_BilingualEduc <- group_by(SV_AllRules, Rule, SV1_Winner_BilingualEduc) %>%
  summarize(SV1_Winner_BilingualEduc_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV1_Winner_BilingualEduc == "Majority") %>% mutate(SV1_Winner_BilingualEduc_Freq = 100 - SV1_Winner_BilingualEduc_Freq)

temp1 <- full_join(Summary_Immigration, Summary_BilingualEduc, by = "Rule")
temp2 <- full_join(Summary_TeacherTenure, Summary_PublicBonds, by = "Rule")

SV1_MinorityVictoryRules <- full_join(temp1, temp2, by = "Rule") %>% 
  select(Rule, matches("Freq"))

stargazer(as.data.frame(SV1_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)

# SV2
Summary_Immigration <- group_by(SV_AllRules, Rule, SV2_Winner_Immigration) %>%
  summarize(SV2_Winner_Immigration_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV2_Winner_Immigration == "Majority") %>% mutate(SV2_Winner_Immigration_Freq = 100 - SV2_Winner_Immigration_Freq)

Summary_PublicBonds <- group_by(SV_AllRules, Rule, SV2_Winner_PublicBonds) %>%
  summarize(SV2_Winner_PublicBonds_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV2_Winner_PublicBonds == "Majority") %>% mutate(SV2_Winner_PublicBonds_Freq = 100 - SV2_Winner_PublicBonds_Freq)

Summary_TeacherTenure <- group_by(SV_AllRules, Rule, SV2_Winner_TeacherTenure) %>%
  summarize(SV2_Winner_TeacherTenure_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV2_Winner_TeacherTenure == "Majority") %>% mutate(SV2_Winner_TeacherTenure_Freq = 100 - SV2_Winner_TeacherTenure_Freq)

Summary_BilingualEduc <- group_by(SV_AllRules, Rule, SV2_Winner_BilingualEduc) %>%
  summarize(SV2_Winner_BilingualEduc_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(SV2_Winner_BilingualEduc == "Majority") %>% mutate(SV2_Winner_BilingualEduc_Freq = 100 - SV2_Winner_BilingualEduc_Freq)

temp1 <- full_join(Summary_Immigration, Summary_BilingualEduc, by = "Rule")
temp2 <- full_join(Summary_TeacherTenure, Summary_PublicBonds, by = "Rule")

SV2_MinorityVictoryRules <- full_join(temp1, temp2, by = "Rule") %>% 
  select(Rule, matches("Freq"))

stargazer(as.data.frame(SV2_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)

# Print them out to the folder 
sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Percent of Minority Victories (per Issue).txt")
stargazer(as.data.frame(SV1_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(SV2_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()



######################################################################
####### Freq Table of Efficient Minority Wins (byRule & Issue) #######
#####################################################################%

# We want two things:
# - proportion of times in which a minority win would have been efficient
# - proportion of times in which a minority win would have been efficient AND was achieved by BV

#--- Prop of Times a Minority Win is Efficient

#--- Prop of Times a Minority Win is Efficient
temp_BE <- SV_RuleAResults %>% group_by(MinorityWinIsEff_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_BE)
temp_IM <- SV_RuleAResults %>% group_by(MinorityWinIsEff_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_IM)
temp_PB <- SV_RuleAResults %>% group_by(MinorityWinIsEff_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_PB)
temp_TT <- SV_RuleAResults %>% group_by(MinorityWinIsEff_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_TT)

MinorityWinEff <- SV_RuleAResults %>% group_by(MinorityWinIsEff) %>% summarise(Count = n()) %>% mutate(Prop = Count/sum(Count)) %>%
  full_join(temp_BE, by = "MinorityWinIsEff") %>% full_join(temp_IM, by = "MinorityWinIsEff") %>%
  full_join(temp_PB, by = "MinorityWinIsEff") %>% full_join(temp_TT, by = "MinorityWinIsEff")
MinorityWinEff[is.na(MinorityWinEff)] <- 0


#--- Prop of Times a Minority Win is achieved by SV

# SV1
temp_BE <- SV_AllRules %>% group_by(Rule, SV1_MinorityWon_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinWithSV = SV1_MinorityWon_BE)
temp_IM <- SV_AllRules %>% group_by(Rule, SV1_MinorityWon_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinWithSV = SV1_MinorityWon_IM)
temp_PB <- SV_AllRules %>% group_by(Rule, SV1_MinorityWon_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinWithSV = SV1_MinorityWon_PB)
temp_TT <- SV_AllRules %>% group_by(Rule, SV1_MinorityWon_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinWithSV = SV1_MinorityWon_TT)

MinorityWinBySV1 <- SV_AllRules %>% group_by(Rule, SV1_MinorityWon) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "SV1") %>% rename(MinorityWinWithSV = SV1_MinorityWon) %>%
  full_join(temp_BE, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule"))

# SV2

temp_BE <- SV_AllRules %>% group_by(Rule, SV2_MinorityWon_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinWithSV = SV2_MinorityWon_BE)
temp_IM <- SV_AllRules %>% group_by(Rule, SV2_MinorityWon_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinWithSV = SV2_MinorityWon_IM)
temp_PB <- SV_AllRules %>% group_by(Rule, SV2_MinorityWon_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinWithSV = SV2_MinorityWon_PB)
temp_TT <- SV_AllRules %>% group_by(Rule, SV2_MinorityWon_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinWithSV = SV2_MinorityWon_TT)

MinorityWinBySV2 <- SV_AllRules %>% group_by(Rule, SV2_MinorityWon) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "SV2") %>% rename(MinorityWinWithSV = SV2_MinorityWon) %>%
  full_join(temp_BE, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("MinorityWinWithSV" = "MinorityWinWithSV", "Rule" = "Rule"))

# Merged 

MinorityWinBySV <- bind_rows(MinorityWinBySV1, MinorityWinBySV2) %>% filter(MinorityWinWithSV == T)
MinorityWinBySV[is.na(MinorityWinBySV)] <- 0

#--- Prop of Times a Minority Win is Efficient AND was achieved by SV

# SV1

temp_BE <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV1_BE) %>% summarise(BE_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV1_BE)
temp_IM <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV1_IM) %>% summarise(IM_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV1_IM)
temp_PB <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV1_PB) %>% summarise(PB_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV1_PB)
temp_TT <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV1_TT) %>% summarise(TT_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV1_TT)

EffMinorityWinBySV1 <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV1) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "SV1") %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV1)  %>%
  full_join(temp_BE, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule"))

# SV2

temp_BE <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV2_BE) %>% summarise(BE_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV2_BE)
temp_IM <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV2_IM) %>% summarise(IM_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV2_IM)
temp_PB <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV2_PB) %>% summarise(PB_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV2_PB)
temp_TT <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV2_TT) %>% summarise(TT_Count = n()) %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV2_TT)

EffMinorityWinBySV2 <- SV_AllRules %>% group_by(Rule, EffMinorityWinWithSV2) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "SV2") %>% rename(EffMinorityWinWithSV = EffMinorityWinWithSV2)  %>%
  full_join(temp_BE, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("EffMinorityWinWithSV" = "EffMinorityWinWithSV", "Rule" = "Rule"))
# Merged 

EffMinorityWinBySV <- bind_rows(EffMinorityWinBySV1, EffMinorityWinBySV2) %>% filter(EffMinorityWinWithSV == T)
EffMinorityWinBySV[is.na(EffMinorityWinBySV)] <- 0


# Print them out to the folder
sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Freq of Efficient Minority Victories.txt")
stargazer(as.data.frame(MinorityWinEff), title = "Prop of times in which Min win is efficient", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(MinorityWinBySV), title = "Prop of times in which SV leads the minority to win (need not be eff)", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(EffMinorityWinBySV), title = "Prop of times in which Min win is efficient and was achieved with SV", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()


######################################
####### Potential Welfare Gain #######
#####################################%

bins25 <- data.frame(Order = seq(0,25,1))
bins25 <- bins25 %>% mutate(Bins = cut(Order, seq(0,25,1), include.lowest = F))

# Potential Welfare Gain: it is the same for all rules, only rule A is used
# Defined as [W(eff) - w(maj)]/w(maj) for all sims
PotentialWelfareGain <- SV_RuleAResults %>%  
  mutate(WelfareFrontier = (AggregateWelfare_Eff - AggregateWelfare_noBV)/AggregateWelfare_noBV) %>%
  select(Round, AggregateWelfare_Eff, AggregateWelfare_noBV, WelfareFrontier) %>%
  mutate(WelfareFrontier = round(WelfareFrontier*100, digits = 2)) %>%
  mutate(Bins = cut(WelfareFrontier, seq(0,25,1), include.lowest = F)) %>%
  group_by(Bins) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) 

UnbinnedPotentialWelfareGain <- SV_RuleAResults %>%  
  mutate(WelfareFrontier = (AggregateWelfare_Eff - AggregateWelfare_noBV)/AggregateWelfare_noBV) %>%
  select(Round, AggregateWelfare_Eff, AggregateWelfare_noBV, WelfareFrontier) %>%
  mutate(WelfareFrontier = round(WelfareFrontier*100, digits = 2)) %>%
  mutate(Bins = cut(WelfareFrontier, seq(0,25,1), include.lowest = F)) 

l1 <- round(mean(na.omit(UnbinnedPotentialWelfareGain)$WelfareFrontier), digits = 2)
l2 <- round(median(na.omit(UnbinnedPotentialWelfareGain)$WelfareFrontier), digits = 2)
l1 <- paste("Mean PWG (% Maj. Welfare): ", l1, sep = "")
l2 <- paste("Median PWG (% Maj. Welfare): ", l2, sep = "")

temp <- left_join(bins25, PotentialWelfareGain, by = "Bins") %>% filter(Order != 0)
temp[is.na(temp)] <- 0


welfarePlot <- ggplot(data = temp, aes(x=Bins, y = Freq)) + geom_bar(stat = "identity") +
  ggtitle("Relative Freq. of Potential Welfare Gain as Percentage of Simple Majority Welfare (No Recal)") + 
  xlab("Potential Welfare Gain (as % of maj. welfare)") + ylab("Relative Frequency") +
  theme(axis.text.x=element_text(size=12, vjust = -.05, angle = 270), title=element_text(size=12)) + 
  annotate("text", x = 18, y = c(.1,.09), label = c(l1,l2)) + coord_cartesian(ylim = c(0, .12)) + 
  scale_x_discrete(drop=FALSE)
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Potential Welfare Gain as Fraction of Simple Majority Welfare (No Recal).png", height = 5.5, width = 9)




######################################
####### Realized Welfare Gains #######
#####################################%

# (a) Majority: W(M)/W(eff) for each simulation. And average over all simulations.
# (b) SV: W(SV)/W(eff) for each simulation. Then both average over all simulations, and average over all s(m, BV) only.

#--- Modified ---#

getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}

NSV_RealizedWelfareGain <- SV_RuleAResults %>%
  mutate(RealWelfGain_All = AggregateWelfare_noBV/AggregateWelfare_Eff) %>%
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5)) %>%
  mutate(Vote = "Simple Maj.")

# NSV_RealizedWelfareGain <- SV_RuleAResults %>%
#   mutate(RealWelfGain_All = AggregateWelfare_noBV/AggregateWelfare_Eff) %>%
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All)) %>%
#   mutate(Vote = "Simple Maj.")

SV1_RealizedWelfareGain <- SV_AllRules %>%
  mutate(RealWelfGain_All = AggregateWelfare_SV1/AggregateWelfare_Eff,
         RealWelfGain_MinSV = ifelse(SV1_MinorityWon == T, AggregateWelfare_SV1/AggregateWelfare_Eff, NA))  %>%
  group_by(Rule) %>% 
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5),
            MeanRealWelfGain_MinSV = mean(RealWelfGain_MinSV, na.rm = T), 
            MinSV_lower95CI = getNthEmpPercentile(RealWelfGain_MinSV,2.5),
            MinSV_upper95CI = getNthEmpPercentile(RealWelfGain_MinSV,97.5)) %>%
  mutate(Vote = "SV1")

# SV1_RealizedWelfareGain <- SV_AllRules %>%
#   mutate(RealWelfGain_All = AggregateWelfare_SV1/AggregateWelfare_Eff,
#          RealWelfGain_MinSV = ifelse(SV1_MinorityWon == T, AggregateWelfare_SV1/AggregateWelfare_Eff, NA))  %>%
#   group_by(Rule) %>% 
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All),
#             MeanRealWelfGain_MinSV = mean(RealWelfGain_MinSV, na.rm = T)) %>%
#   mutate(Vote = "SV1")

SV2_RealizedWelfareGain <- SV_AllRules %>%
  mutate(RealWelfGain_All = AggregateWelfare_SV2/AggregateWelfare_Eff,
         RealWelfGain_MinSV = ifelse(SV2_MinorityWon == T, AggregateWelfare_SV2/AggregateWelfare_Eff, NA))  %>%
  group_by(Rule) %>% 
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5),
            MeanRealWelfGain_MinSV = mean(RealWelfGain_MinSV, na.rm = T), 
            MinSV_lower95CI = getNthEmpPercentile(RealWelfGain_MinSV,2.5),
            MinSV_upper95CI = getNthEmpPercentile(RealWelfGain_MinSV,97.5)) %>%
  mutate(Vote = "SV2")

# SV2_RealizedWelfareGain <- SV_AllRules %>%
#   mutate(RealWelfGain_All = AggregateWelfare_SV2/AggregateWelfare_Eff,
#          RealWelfGain_MinSV = ifelse(SV2_MinorityWon == T, AggregateWelfare_SV2/AggregateWelfare_Eff, NA))  %>%
#   group_by(Rule) %>% 
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All),
#             MeanRealWelfGain_MinSV = mean(RealWelfGain_MinSV, na.rm = T)) %>%
#   mutate(Vote = "SV2")


AvgRealizedWelfareGain <- bind_rows(SV1_RealizedWelfareGain, SV2_RealizedWelfareGain, NSV_RealizedWelfareGain)


sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Average Realized Welfare Gain (per Rule).txt")
stargazer(as.data.frame(AvgRealizedWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
sink()

AvgForegoneWelfareGain <- AvgRealizedWelfareGain %>% 
  mutate(MeanForegoneGain_All = 1 - MeanRealWelfGain_All, MeanForegoneGain_MinSV = 1 - MeanRealWelfGain_MinSV) %>%
  select(-matches("Real")) %>% select(-matches("CI"))

# AvgForegoneWelfareGain <- AvgRealizedWelfareGain %>% 
#   mutate(MeanForegoneGain_All = 1 - MeanRealWelfGain_All, MeanForegoneGain_MinSV = 1 - MeanRealWelfGain_MinSV) %>%
#   select(-matches("Real"))

sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Average Foregone Gain (per Rule).txt")
stargazer(as.data.frame(AvgForegoneWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
sink()


#--- Modified (end) ---#

###################################################
####### Freq. of Welfare Increases (byRule) #######
##################################################%

# SV1 - (Fig. 5 upper half; in Elections paper)

SV1_AllRules <- filter(SV_AllRules, SV1_MinorityWon == TRUE) %>% 
  mutate(SV1_WelfareGain = ifelse(AggregateWelfare_SV1 > AggregateWelfare_noBV, TRUE, FALSE)) %>%
  group_by(Rule, SV1_WelfareGain) %>% summarize(Count = n()) %>% 
  mutate(Freq = round(Count/sum(Count), digits = 3)*100) %>% filter(SV1_WelfareGain == TRUE)

welfarePlot <- ggplot(data = SV1_AllRules, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("SV1 (No Recal) - Freq. of Welfare Increases") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Welfare Increases") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV1 (No Recal) - Freq. of Welfare Increases.png", height = 5.5, width = 5.5)

#SV2 - (Fig. 5 upper half; in Elections paper)

SV2_AllRules <- filter(SV_AllRules, SV2_MinorityWon == TRUE) %>%
  mutate(SV2_WelfareGain = ifelse(AggregateWelfare_SV2 > AggregateWelfare_noBV, TRUE, FALSE)) %>%
  group_by(Rule, SV2_WelfareGain) %>% summarize(Count = n()) %>% 
  mutate(Freq = round(Count/sum(Count), digits = 3)*100) %>% filter(SV2_WelfareGain == TRUE)

welfarePlot <- ggplot(data = SV2_AllRules, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("SV2 (No Recal) - Freq. of Welfare Increases") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Welfare Increases") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV2 (No Recal) - Freq. of Welfare Increases.png", height = 5.5, width = 5.5)

# table of results (combined)

FreqWelfareIncrease <- bind_rows(mutate(SV1_AllRules, Vote = "SV1"), mutate(SV2_AllRules, Vote = "SV2")) %>% select(-matches("Welfare"))
sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Freq. of Welfare Gains (per Rule).txt")
stargazer(as.data.frame(FreqWelfareIncrease), type = "text", summary = FALSE, rownames = FALSE, digits = 1)
sink()


###################################################
####### Mean % Welfare Change Plot (byRule) #######
##################################################%

# SV1 -  (Fig. 5 bottom half; in Elections paper)
SV1_AllRules <- filter(SV_AllRules, SV1_MinorityWon == TRUE) %>% 
  mutate(SV1_PercentWelfareChange = (AggregateWelfare_SV1 - AggregateWelfare_noBV)/AggregateWelfare_noBV) %>%
  group_by(Rule) %>% summarize(SV1_MeanPercentWelfareChange = round(mean(SV1_PercentWelfareChange), digits = 3)*100) %>%
  rename(MeanPercentWelfareChange = SV1_MeanPercentWelfareChange)

welfarePlot <- ggplot(data = SV1_AllRules, aes(x=Rule, y = MeanPercentWelfareChange)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5)) +  
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Welfare Change") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV1 (No Recal) - Bootstrapped Mean Percent Welfare Change")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV1 (No Recal) - Mean Percent Welfare Change.png", height = 5.5, width = 5.5)


# SV2 -  (Fig. 5 bottom half; in Elections paper)
SV2_AllRules <- filter(SV_AllRules, SV2_MinorityWon == TRUE) %>% 
  mutate(SV2_PercentWelfareChange = (AggregateWelfare_SV2 - AggregateWelfare_noBV)/AggregateWelfare_noBV) %>%
  group_by(Rule) %>% summarize(SV2_MeanPercentWelfareChange = round(mean(SV2_PercentWelfareChange), digits = 3)*100) %>%
  rename(MeanPercentWelfareChange = SV2_MeanPercentWelfareChange)

welfarePlot <- ggplot(data = SV2_AllRules, aes(x=Rule, y = MeanPercentWelfareChange)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5)) +  
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Welfare Change") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV2 (No Recal) - Bootstrapped Mean Percent Welfare Change")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/SV2 (No Recal) - Mean Percent Welfare Change.png", height = 5.5, width = 5.5)

# table of results (combined)

MeanWelfareChange <- bind_rows(mutate(SV1_AllRules, Vote = "SV1"), mutate(SV2_AllRules, Vote = "SV2")) 
sink("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Table - Mean Percent Welfare Change (per Rule).txt")
stargazer(as.data.frame(MeanWelfareChange), type = "text", summary = FALSE, rownames = FALSE, digits = 2)
sink()


###################################################
####### Hist of % Welfare Change (per Rule) #######
##################################################%

WelfareChangeHist <- function(sampleData, SetRule, SampleTitle)
{
  #--- Set Graphing Parameters
  
  bins_width5 <- c("[-20,-15)","[-15,-10)","[-10,-5)", "[-5,0)", "[0,5)", "[5,10)", "[10,15)", "[15,20)", "[20,25)", "[25,30)")
  bins5 <- data.frame(Bins = bins_width5, BinOrder = seq(1,length(bins_width5),1), stringsAsFactors = F)
  
  #--- Determines the current vote
  if(SampleTitle== "SV1"){
    tempRule <- filter(sampleData, SV1_MinorityWon == TRUE, Rule == SetRule) %>%
      mutate(SV1_PercentWelfareChange = (AggregateWelfare_SV1 - AggregateWelfare_noBV)/AggregateWelfare_noBV,
             SV1_PercentWelfareChange = round(SV1_PercentWelfareChange, digits = 3)*100,
             Bins = cut(SV1_PercentWelfareChange, right = F, breaks = seq(-20, 30, 5)))
  }
  
  if(SampleTitle == "SV2"){
    tempRule <- filter(sampleData, SV2_MinorityWon == TRUE, Rule == SetRule) %>%
      mutate(PercentWelfareChange = (AggregateWelfare_SV2 - AggregateWelfare_noBV)/AggregateWelfare_noBV,
             PercentWelfareChange = round(PercentWelfareChange, digits = 3)*100,
             Bins = cut(PercentWelfareChange, right = F, breaks = seq(-20, 30, 5))) 
  }
  
  tempRule <- tempRule %>% group_by(Bins) %>% summarize(Count = n()) %>% mutate(Freq = round(Count/sum(Count), digits = 3)*100) 
  tempRule <- full_join(bins5, tempRule, by = "Bins")
  tempRule[is.na(tempRule)] <- 0
  tempRule$Bins <- factor(tempRule$Bins, levels = tempRule$Bins[order(tempRule$BinOrder)])
  
  welfarePlot <- ggplot(tempRule, aes(x = Bins, y = Freq)) +
    geom_bar(stat="identity") + xlab("Percent Welfare Change") + ylab("Relative Frequency") +
    scale_y_continuous(limits = c(0,70), breaks = seq(0, 70, by = 10)) +
    theme(axis.text.x=element_text(size=12, vjust = .5, angle = 300), legend.text=element_text(size=16), title=element_text(size=12)) +
    ggtitle(paste("Rule ",SetRule, " (", SampleTitle, ") ", " - Histogram of % Welfare Change (No Recal)", sep = ""))
  print(welfarePlot)
  ggsave(welfarePlot, file=paste("./SV & QV Simulations (No Recalibration)/SV - Welfare & Minority Victories/Rule ", SetRule, " (", SampleTitle,") - Histogram of Percent Welfare Change (No Recal).png",sep=""), height = 4.5, width = 5.5)
}

# SV1 - (Fig. 6 in Elections paper)

WelfareChangeHist(SV_AllRules, "A", "SV1")
WelfareChangeHist(SV_AllRules, "B", "SV1")
WelfareChangeHist(SV_AllRules, "C", "SV1")
WelfareChangeHist(SV_AllRules, "D", "SV1")
WelfareChangeHist(SV_AllRules, "E", "SV1")

# SV2 - (Fig. 6 in Elections paper)


WelfareChangeHist(SV_AllRules, "A", "SV2")
WelfareChangeHist(SV_AllRules, "B", "SV2")
WelfareChangeHist(SV_AllRules, "C", "SV2")
WelfareChangeHist(SV_AllRules, "D", "SV2")
WelfareChangeHist(SV_AllRules, "E", "SV2")




###########################################
##### Aggregate Preference Histograms #####
###########################################


ProposalHist <- function(sampleData, Proposal, RawPrefSummary)
{
  #--- Set Graphing Parameters
  
  plotMax <- 70

  bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                    "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                    "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")
  
  
  # Calculate the average total points for each of the proposals and store it on the dataframe
  PrefSummary <- summarise_all(RawPrefSummary, funs(mean))
  
  # Calculate the average number of subjects and store it on a temp dataframe
  temp <- summarise_all(sampleData, funs(mean))
  temp <- unname(unlist(temp[1,]))
  
  # Create dataframe with summarized aggregate data
  currentProp <- data.frame(Bins = factor(bins_width10), NumSubjects = temp, BinOrder = seq(1,21,1), 
                        BinColor = c(rep(0,10),1,rep(0,10)), stringsAsFactors = F)
  currentProp$Bins <- factor(currentProp$Bins, levels = currentProp$Bins[order(currentProp$BinOrder)])
  
  #plot
  intensityPlot <- ggplot(data = currentProp, aes(x=Bins, y = NumSubjects, fill=BinColor)) + 
    geom_bar(stat="identity") + scale_y_continuous(limits = c(0,plotMax), breaks = seq(0, plotMax, by = 10)) +  
    xlab("Number of Points Assigned to Proposal") + ylab("Number of Subjects") + 
    theme(axis.text.x=element_text(size=12, vjust = -.05, angle = 90), legend.text=element_text(size=16), title=element_text(size=14)) +
    ggtitle(paste("SV (No Recal) - ",Proposal, sep = "")) + guides(fill = F)
  
  l1 <- paste("Mean Total Points InFavor: ",PrefSummary[[paste(Proposal,"_PointsInFavor", sep = "")]], sep = "")
  l2 <- paste("Mean Total Points Against: ",PrefSummary[[paste(Proposal,"_PointsAgainst", sep = "")]], sep = "")
  intensityPlot <- intensityPlot + annotate("text", x = 18, y = c(70,66), label = c(l1,l2))
  
  print(intensityPlot)
  ggsave(intensityPlot, file=paste("SV & QV Simulations (No Recalibration)/SV - Preference Histograms/SV (No Recal) - ",Proposal,".png", sep = ""), height = 4.5, width = 9.5)
  
}

ProposalHist(sampleData = SV_ImmBSPref, Proposal = "Immigration", RawPrefSummary = SV_PrefSummary_raw)
ProposalHist(sampleData = SV_BondBSPref, Proposal = "PublicBonds", RawPrefSummary = SV_PrefSummary_raw)
ProposalHist(sampleData = SV_EducBSPref, Proposal = "BilingualEduc", RawPrefSummary = SV_PrefSummary_raw)
ProposalHist(sampleData = SV_TeachBSPref, Proposal = "TeacherTenure", RawPrefSummary = SV_PrefSummary_raw)

